%CD 'U:/'
job 'test programme to obtain maximum force value all parameters to determine modulus of elasticity and modulus of rupture 10_11_10 CLM '

"  CALCULATION of FLEXURAL STIFFNESS (EI)
	EI = (F/V * (a^2/12)*(3L-4a)

	where F/V = maximumum slope (Smax)
	a = distance from load point to support (35mm)
	L = length between 2 supported ends (70 mm)"

"  CALCULATION ADJUSTMENT for CROSS SECTIONAL AREA
	I=(3.142*(D2**4-D1a**4))/64
	where
	D2=(2*(sqrt(stem area/3.142)))
	D1a=(2*(sqrt(hollow area area/3.142)))"
	
	




"  CALCULATION of MODULUS of ELASTICITY
	EI/I"

"  CALCULATION of MODULUS of RUPTURE
	MOR=(Fmax*a*(D2/2-D1a/2))/I				 
	where Fmax= peak force
	a=distance from load point to support (35mm)
	D2=(2*(sqrt(stem area/3.142)))
	D1a=(2*(sqrt(hollow area area/3.142)))
	I(3.142*(D2**4-D1a**4))/64"
	
	





output [width=170]1		"sets outout width of screen (=1)"

text filename,Filename
scalar t1,t2,slope,seslope,tslope,Fmax,Smax
calc Smax=0
" insert locations of input files and desired location for calculation outputs"

open 'stress-strail file location.txt';channel=2			"sets input channel for test data to 2"
open 'stem cross-sectional measurements file location.txt';channel=3			"sets input channel for stem data to 3"
open 'output file location.txt';channel=3;filetype=o;width=200	"sets output channel"


print [ch=3;iprint=*]'Test name                 Peak force          time(s)  \
 Max. Slope   at  seconds	Flexural stiffness D2   D1a  Ia MOE MOR'	"headers in output file"

for test=1...205	"sets number loops to equal no of samples"

	read [ch=2;print=*]filename						"reads header data from file 2"
	read [ch=2;print=*]kg,sec
										"reads test data from file 2"
	read [ch=3;print=*]filename,stemarea,hollowarea	
		"reads stem data from file 3"

"dgraph kg;sec"						"outputs to high quality graph of test data"

"	CALCULATES PEAK FORCE (Fmax)	"

restrict kg;sec.lt.15					"limits the test data to first part of output ie < 3 sec"

	sort [index=kg;dir=d]kg,sec;Skg,Ssec		"sorts force into descending order & time indexed"
	calc TIMEfmax=elements(Ssec;1)			"time in sec of Fmax"
	  &  Fmax=elements(Skg;1)			"Peak force (Fmax)"

	calc inc=TIMEfmax/.005				"calculate number of increments to reach Fmax"

restrict kg						"removes limit"

"	CALCULATES SLOPE (Smax (F/V)) and other PARAMETERS (adjust depending on sample)	"

	calc t1=-0.005					"sets 'span' of regression slope in seconds"
	  &  t2=0.025					"time increments on 0.005 sec intervals"
	  &  a=6				"sets load distance"
	  &  L=12
	  &  D2=(2*(sqrt(stemarea/3.142)))
	  &  D1a=(2*(sqrt(hollowarea/3.142)))
	  
	 

	  
	    
				"sets span between supports"

				
	  		

		for i=1...inc				"sets no of increments for regressions"
							
			calc t1=t1+0.005		"sets distance(s)between regressions"
			calc t2=t2+0.005

			restrict kg;(sec.ge.t1).and.(sec.le.t2)	"limits calculation to a regression between time interval = 0.3 s = 6 data points"

				model kg			" Y axis = dependent trait"
				fit [print=*]sec		" X axis = independent trait"
				rkeep kg;meandev=resms;tmeandev=totms;estimates=sl;se=tsl;df=resdf;tdf=totdf	 "stores output from regression"

				equate [oldform=!(-1,1)]sl;slope
				  &    [oldform=!(-1,1)]tsl;seslope

				calc   tslope=slope/seslope
				  &    totSS=totms*totdf				
				  &    resSS=resms*resdf				
				  &    regSS=totSS-resSS
				  &    Rsq=(regSS/totSS)*100

					"	CALCULATION of MAXIMUM SLOPE (Smax)	"

					if slope.gt.Smax
						calc Smax=slope
						  &  time=(t1+t2)/2  			"(*sec)eg 15 * "	
					endif

				"	CALCULATION of MODULUS of ELASTICITY	"

					calc EI=(Smax*(a**2/12)*((3*L)-(4*a)))	"Flexural stiffness"
					  &  I=(3.142*(D2**4-D1a**4))/64	"Second moment of area"
					  &  MOE=EI/I
					
							"Modulus of elasticity"

			restrict kg		"removes limits of regression (t1 & t2)"
					
				"	CALCULATION of MODULUS of RUPTURE	"
		for i=1...inc				
					calc rupture=(Fmax*a*(D2/2-D1a/2))
											"MOR"
  					&    I=(3.142*(D2**4-D1a**4))/64
					&    MOR=(rupture/I)  
					
					
		restrict kg						
							
		endfor

		endfor

	print [ch=3;iprint=*;squash=y]filename,Fmax,TIMEfmax,Smax,time,EI,D2,D1a,I,MOE,MOR;fieldwidth=30,15,15,15,15,15,15,15,15,15,15 "output to file"
	print [iprint=*;squash=y]filename,Fmax,TIMEfmax,Smax,time,EI,D2,D1a,I,MOE,MOR;fieldwidth=30,15,15,15,15,15,15,15,15,15,15	"output to screen"

delete [redefine=yes]sec,kg

calc slope=0		"resets scalars to zero for next loop & set data"
  &  Smax=0
  &  time=0
  &  EI=0
  &  I=0
  & Fmax=0
endfor

close ch=2,3		"close input files after all calculations"

stop




